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We carry out a stability analysis for tiie real space split operator method for the propagation of the time-dependent 
Klein-Gordon equation that has been proposed Ruf et al. [M. Ruf, H. Bauke, C.H. Keitel, A real space split 
operator method for the Klein-Gordon equation, Journal of Computational Physics 228 (24) (2009) 9092-9106, 
"doiilO. 1016/j.jcp.2009.09.012|. The region of algebraic stability is determined analytically by means of a 
von-Neumann stability analysis for systems with homogeneous scalar and vector potentials. Algebraic stability 
implies convergence of the real space split operator method for smooth absolutely integrable initial conditions. 
In the limit of small spatial grid spacings h in each of the d spatial dimensions and small temporal steps t, the 
stability condition becomes h/r > yfdc for second order finite differences and V3/j/(2t) > yfdc for fourth order 
finite differences, respectively, with c denoting the speed of light. Furthermore, we demonstrate numerically that 
the stability region for systems with inhomogeneous potentials coincides almost with the region of algebraic 
stability for homogeneous potentials. 
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1. Introduction 

The Klein-Gordon equation is a partial differential equation that governs the quantum evolution of wave functions for relativistic 
spinless particles. It finds its applications for example in modeling light matter interaction in the relativistic regime |1-3|. 
The derivation of analytical solutions for the time-dependent problems in relativistic quantum dynamics requires sophisticated 
mathematical tools |4-9| and many setups necessitate numerical approaches solving the Dirac equation |10-17| or the Klein- 
Gordon equation 1 16 18|, especially problems with time-dependent potentials. In weakly relativistic regimes the solution of the 
Pauli equation with first order relativistic corrections may be sufficient [19]. 

Ruf et al. II20II developed an explicit real space split operator finite difference scheme for the numerical solution of the time- 
dependent Klein-Gordon equation. It has been applied for example to study relativistic ionization characteristics of laser-driven 
hydrogen like ions 1 16| and numerical signatures of non-self -adjointness in quantum Hamiltonians [21 J. As an explicit scheme, the 
real space split operator finite difference scheme allows for an efficient implementation on distributed memory parallel computers. 
Explicit finite difference schemes for hyperbolic partial differential equations, however, are known to be never unconditionally 
stable 1 22 23). This means, that the numerical eiTor that is caused by the application of an explicit finite difference scheme to 
the Klein-Gordon equation will grow exponentially unless the spatial and temporal grid spacings fulfill certain conditions. The 
purpose of this contribution is to study the stability of the explicit real space split operator finite difference scheme as introduced 
in 1 20 1 and to derive stability conditions for this scheme. 

The reminder of the paper is organized as follows. In order to keep the paper self-contained, we will give a brief description of 
the real space split operator method for the Klein-Gordon equation in Sec.|2j Section[3]wiII state more precisely the notion of 
numerical stability by defining strong and algebraic stability, while Sec.|4]reviews conditions for these notions of stability. A von- 
Neumann stability analysis for constant potentials will be performed in Sec.[5]and in Sec.[6]we will compare our theoretical results 
with the stability of an actual computer implementation of the real space split operator method for systems with homogeneous 
potentials as well as for systems with spatial dependent potentials. 



2. A real space split operator method for tlie Klein-Gordon equation 

The Klein-Gordon equation — sometimes also called Klein-Fock-Gordon or Klein-Gordon-Schrodinger equation — is an equation 
of motion for a complex- valued scalar quantum wave function ip(x, t) Il8l l24ll25l . It governs the behavior of a charged spinless 
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particle with mass m and charge q moving under the effect of electrodynamic potentials A{x, t) and (p{x, t). Introducing the speed 
of light c and the reduced Planck constant h, the scalar Klein-Gordon equation reads 



d_ 



ih^ - q(p{x, t)\ -c^ (-i/iV - qA{x, t)f - m^c'^ 



(f(x,t)^0. (1) 



This hyperbolic partial differential equation is of second order in time. For numerical purposes a transformation of ([T]) into a first 
order partial differential equation is advantageous. Introducing the two-component wave function L24i 



-''Mitt 



( { {^{X, t)+^ {inf, - q4>(x, t)) ipix, t)) 
\ [^{X, t)-^ (iS| - q<p(x, t)) ifix, t)l 



(2) 



the Klein-Gordon equation ([T]) may be transformed into Schrodinger form with the Hamiltonian operator H, viz. 

^j^ d'i'Oct) ^ fj^jy^^^^ ^ lo^+Wji _ ^^^^^ ^^^2 ^ ^^(^^ ^ o-^niA f) . (3) 

at \ 2m j 

The wave-function's components are related by the complex Hermitian unitary Pauli matrices 0-2 and ctt, that obey the Pauli 
algebra 

o-iO-j = ieijMCTk + Si J , (4a) 
[cr;,cry] = 2ie,j;tcrj. , (4b) 
{cri,CTj}^26ij, (4c) 

with i, j,k e {1,2,3} and where eij^k denotes the permutation symbol and 6ij is the Kronecker delta. The Pauli algebra does not 
define the Pauli matrices uniquely. In numerical applications, we choose the representation 

^^' = (1 0)' '^^^(i 0)' '^^^"(o -i)- 

Solutions ^{x, t) of equation ^ satisfy the continuity equation 

dp{x, t) 

-^^+V-y(jc,f) = (6) 

with the density 

p{x, t) = "V^ix, t)(Ti'V{x, t) = \'i',{x, tt - Wlix, tf . (7) 

and the current 

j{x, t) = [T*(jc, t)(T^{(T^ + icr2)V>I'(jc, t) - {V^\x, f))cr3(cr3 + i(T2mx, t)] - '^^^^*{x, t)a^{(T2 + i(T2mx, t) , (8) 
Im m 

where *I'*(jc, t) denotes the complex conjugate of *P(j!:, f)- The density (|7]i is not positive definite and, therefore, the integral 

\mx,ml^^ J^'*(x,t)cr3'i'(x,t)d^ (9) 

over the t/-dimensional space is not a norm and p(x, t) cannot be interpreted as a probability density and j{x, t) is not a probability 
density current. The quantity qp{x, t), however, may be interpreted as a charge density and q j(x, t) as a charge current Il5l48l l24ll . 
The conservation of the pseudo norm Q implies charge conservation. Note that the Euclidean norm 

,im«ii= = /.'-<.,omod^ (,o, 



is not conserved under the evolution of the Klein-Gordon equation. 

In order to derive the real space split operator method of |l20l, we start from the formal solution of (|3]l with the initial condition 
^{x, 0) that may be expressed as 

'V{x,t)^U{t,QY¥{x,Q) (11) 
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by introducing the unitary time-evolution operator 

U{t2,h) = f exp 



(12) 



and Dyson's time ordering operator t. The time-evolution operator U{t2, t\) effects the wave-function's evolution from time t\ to 
time t2- Let 0(t) denote some possibly time-dependent operator and define the operator 



UdihJuS) = exp J-(5^ J^' d(f')df'j , 



(13) 



that depends on the times t\ and t2 and the auxiliary parameter 6. Expanding U (t + t, t) to the third order in t and introducing the 
operators 



era H- icT-, , 

K{x, t) = ^- i-im - qA{x, t)f 

Zm 

V{x, t) - qip{x, t) + cr^mc^ , 



(14a) 
(14b) 



the time-evolution operator for the Klein-Gordon equation ( 12 1 can be factorized as 



= Uy[t + T,t,\)U^(t + T,t,l)Uy[t + T,t,\) + 0[T^) . 



(15) 



Neglecting terms of order O (t^), equation jisj gives an explicit second order accurate time-stepping scheme for the propagation 
of the Klein-Gordon wave function ^'(x, t). 

The operator Uy{t + T,t, 6) is diagonal in real space and explicitly given by 



Uy it + T,t,6)^'(X,t) = 



exp (-^i j''^^ q(f>(x, t') + mc^ df') *Pi (x, f)' 
exp J^'^^ q^{x, f) - mc^ df) 'V2(x, t), 
The action of the operator U^it + T,t, 5) on ^{x, t) may be calculated by the Taylor series of the exponential function 

■HI 



(16) 



it + T, t, 6) 'i'ix, t) = exp 



K{x,t')dt'\^>(x,f) 



(17) 



Recognizing that all but the first two terms of the infinite series ( 17 1 vanish because the operator (1x3 + 10-2) is nilpotent, that is, 
(0-3 + \0-2f' - 0, yields 



Uk it + T, t, 6) "fix, t) = "Vix, t) - (0-3 + xa-2) 



6\ 
2mh 



V-qA{x,t')\ '¥(x,t)dt' 



Introducing the quantities 



Co = exp 



2n j, 



<f>(x, t') df' , 



c\ 

C2 

C3 



2m ' 

- f A(x,t')dt' , 
m J, 



— V ■ A(x, f) - —A(x, t'f df' . 
2m 2mh 



(18) 

(19a) 
(19b) 
(19c) 

(19d) 
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and taking advantage of the standard representation ^ of the Pauh matrices the operator product ( 15 1 finally may be transformed 
into a product of the compact matrices 



U(t + T,t) = 



(*) 

Co 0\(e 
Co 



-imc^T ji2h) 





^-\-imc-Tj{2ti) 



1 + (CiV^ + C2 ■ V + C3) + (ciV^ + C2 ■ V + C3) 
- (CiV2 + C2 • V + C3) 1 - (CiV2 + C2 ■ V + C3) 

Co 0\/e-™'''^/(2fi) 



CO 









-^+imc- T I (Ih) 



0[t'). (20) 



In a computer implementation of the real space split operator method, the wave function T(jc, f) is discretized on a rectangular 
id + l)-dimensional space-time lattice with spatial grid spacings hj (i e {1, . . . , d}) and time steps of size t, viz. 



hitii 



hdnd. 



with e 



f - Tk with k e 



(21a) 



(21b) 



The quantity h will refer to a t/-dimensional vector with components /i,. The discretized two-component wave function will be 
denoted by 



2,n 



XTJK 



(22) 



and identifies a vector with components T„ for fixed k. The discrete version of the pseudo norm (|9|) reads 

d 



(23) 



1=1 



For the stability analysis we will also need the proper norm 



(24) 



!=1 



First and second order derivatives in ([20)i are approximated by finite differences. We will utilize the second order approximations 



V2,„/(x) and ,/(x) 



dfjx) 
dx 



V2,/,/(x) + (9(/i2) 



-f(x - h) + f(x + h) 
2h 



5x2 

as well as the fourth order schemes V4_/,/(x) and i^f{x) 



5/(x) 



V4,/,/(x) + 0(/i4) 



fix - 2h) - 8/(x - h) + 8/(x H- /;) - f{x + 2h) 



dx ' ^ ^ 12h 

d^fix) „2 , , ^ /, 4\ -fix - 2h) + I6f(x -h)- 30f(x) + 16/(x + h) - fix + 2h) 



dx^ 



= V2,/(x) + 0(//) = 



12/!2 



(25a) 
(25b) 

(26a) 
(26b) 



Plugging (25 I or (26 1 into (20 1 and taking into account boundary conditions yields the finite difference scheme (the discretized 
version of {T5]l) 



(27) 



with some sparse 2N x 2N matrix Uh,T that depends on the grid spacings h and r, where denotes the total number of spatial 
grid points. Thus, the real space split operator method for the Klein-Gordon equation is essentially a matrix vector multiplication. 
A stability analysis has to examine the properties of this matrix vector multiplication as a function of h and t. Before we can 
proceed with a stability analysis, however, we have to give the term numerical stability a precise definition. 
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3. Strong and algebraic stability 



Loosely speaking numerical stability of a finite difference scheme means that numerical errors that are caused by the discretization 
of a partial differential equation are bounded. There is no unique way to cast this intuitive notion into a formal mathematical 
statement and different mathematical concepts of stability may be found in the literature. We will consider the so-called strong 
stability Il27ll28i and algebraic stability Il28ll29l . An explicit finite difference scheme ^f'^^' - Uh,T^'' for a partial differential 
equation, that is first-order in time, is called strongly stable in the stability region A if there is for any fixed positive time t a 
constant C such that 



hp 



for all (h,T)eA, < kr < t . 



(28) 



Strong stability requires that the norm is bounded independently of the grid spacings t and h. An explicit finite difiference 

scheme *P*^' = Uhj^'' for a first-order partial differential equation is called algebraically stable in a stability region A if there are 
for any fixed positive time f two constants C and a > such that 



< Ct"" ||*P°||^ for all (h,T)eA, < kr < t . 



(29) 



Algebraic stability allows the norm to grow proportionally to r"" and, therefore, it is not necessarily bounded in the limit r — > 0. 
Note that algebraic stability is a rather weak stability concept; there is no intermediate regime between algebraic stability and 
exponential explosion. Strong stability implies algebraic stability. For partial differential equations for scalar unknown functions 
strong and algebraic stability are equivalent. For systems of partial differential equations as the Klein-Gordon equation 
however, strong stability does not follow by algebraic stability ESl . 

Strong and algebraic stability are closely related to convergence. A finite difference scheme is called convergent if the 
discretization error satisfies 



\\'i'(x,Tk)-'i'\^0 for h,T^O. 



(30) 



Under certain conditions, strong and algebraic stability imply convergence. A consistent finite difference scheme for a partial 
differential equation for which the initial value problem depends continuously on the initial condition and that is of first order in 
its time derivative is convergent if and only if it is strongly stable |27 |. An explicit finite difference scheme = Uh j^'^ for a 
first order partial differential equation '\Tid^(x, t)/dt = H{t)^'{x, t) is called consistent if 



>P*+' - »P* Uh.T - 1 







for h,T^O. 



(31) 



Note that a consistent strongly stable scheme is convergent independently of the initial condition. If a finite difference scheme 
4^*^+' = t//l,T*^'* is algebraically stable and Uhj does not depend on spatial coordinates, then it is convergent provided that the 
Fourier transform *P(/;, 0) of the initial condition "^{x, 0) decays faster than an appropriate power \p\^' for \p\ oo [28 J. Note that 
this condition applies only to homogeneous problems (Uh,T is not allowed to depend on spatial coordinates.) and it is a sufficient 
condition for convergence but not necessary. The condition on the initial function ^{x, 0) is satisfied by any function that is at 
least /-times differentiable and its fth derivative is absolutely integrable. 



4. Stability conditions 



There are various necessary and/or sufficient criteria 



for strong or algebraic stability of finite difference schemes of 



the type ( [27| . For example, one can show that a finite difference scheme of the type \21) is strongly stable in a region A if and 

(32) 



only if for every time /, there exists a constant C such that powers of Uhj are bounded as 



\U 



h,T\ 



< C for all (h,T)eA, < kr < t . 



For a stability analysis in Fourier space one phrases the discrete function ^ via its Fourier transform ^ (p), viz. 



(33) 



Plugging the ansatz ( |33| l into a finite difference scheme of type ( [27| for a partial differential equation with constant coefficients 
and assuming a spatial grid of infinite size individual Fourier modes decouple yielding a propagation equation 



(34) 
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for each Fourier mode ^ e {-n, nY- The equivalent to the condition ( 32 1 in Fourier space reads 

||f7/,,T(^)*|| < C for all ^ e {-n, nY , (h,T) e A, < kr < t . 



(35) 



The condition ( |35| ) follows from ( (32] i via Parseval's Theorem. 

Proving the boundedness of powers of a matrix is often not easy, thus other criteria may be applied. The von-Neumann criterion 
requires for a one-step finite difference scheme with constant coefficients (27 1 that there is a constant K such that for all (h, r) in 
the stability region A the eigenvalues UhjiO of Uhji^) are bounded as 



M/l. 



r{^)\<\+KT. 



(36) 



The von-Neumann stability criterion provides a sufficient condition for algebraic stability f28l. Because strong stability implies 
algebraic stability, the von-Neumann stability criterion is a necessary criterion for strong stability but in general not sufficient. In 
the stability analysis of Sec.|5]we will use the slightly stronger criterion [m^ ,-(^)| ^ 1- 



5. Von-Neumann stability analysis for homogeneous potentials 



In this section we are going to perform a von-Neumann stability analysis for the real space split operator method for the 
Klein-Gordon equation with homogeneous potentials 



(p{x, f) = (f>i) A{x, t) = Ac 



(37) 



Transforming the matrix d20b with finite difi^erences of second order ( 25 1 or fourth order (p6]l into Fourier space the finite difference 



operators in ( 20 1 have to be replaced by 



^Xhfix) ^ V3,ft/(^) = -i 



sin(frf) 

d 



i=l 



2-2 cos(^,) 

^2 



or 



/ 8sin(fi)-sm(2^i) N 
6hi 



!=1 



15- 16cos(^,)-(-cos(2^,) 



respectively. With ( |38) and ( |39] l the Fourier space propagation matrix becomes for homogeneous potentials 

y ^ 2m n.n m 2mn j j \ 2m n,n m 2mn j 



(38a) 



(38b) 



(39a) 



(39b) 



(40) 



where n 6 {2,4). To derive (40i we have used that for a homogeneous scalar potential the matrix (*) in (20 1 commutes with all 
other matrices in 1 20 1. The moduli of the eigenvalues of (40 1 do not depend on the factor Q-^i'PoTlh^ Therefore, for a von-Neumann 
stability analysis we can set = without loss of generality. The numerical stability does not depend on the magnitude of the 
scalar potential. For convenience, we introduce the real valued quantities 



K{^) = y Kii^d with Kii^d = A V2 + ^ V,,A, - 

^—i 7.m ' ' im 



2m 



mi 



2mh 



(41) 
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2 4 6 8 10 2 4 6 8 10 

mch/h mch/ti 

FIG. 1 : Regions of algebraic stability A (shaded in gray) for a one-dimensional finite difference scheme with three point finite differences \25\ 
for two different values of the homogeneous vector potential Aq as a function of the spatial grid spacing h and the temporal step width r. The 
larger the vector potential the smaller the stability region. 



and 



y(^) = cos (ntc^T/fij + k{^)t sin (mc^r/tij 
such that the eigen-values of (|40|t are given by 



(42) 



(43) 



The quantity y(^) is real valued and, therefore, the von-Neumann criterion |m/,,t(^)| < 1 is equivalent to y(^)^ < 1 which must be 
fulfilled for every ^ e [-7t,7t]'^. The von-Neumann criterion is fulfilled for all ^ e [-tt, tt]'' if it is fulfilled for all ^* that make 
yi^)^ extremal. Qualification for extremality of y(^)^ is 



2y(^*)T sin (jnc^r/lij 



for all /€ {!,...,£/). 



(44) 



Thus, to maximize yi^) 



d^i 



forall/e {!,...,£/) 



(45) 



is required. The region A of algebraic stability for the scheme (27 1 with constant potentials is finally given by 



A = \ih,T)e [0, oo)''+i with yiff < 1, for all f e [-n, nf that are the solutions of 



dKii^i) 



d^i 



(46) 



Figure [T] shows the region of algebraic stability A for a one-dimensional finite difference scheme with three point finite 
differences ( [25] l for two different values of the homogeneous vector potential as a function of the spatial grid spacing h and the 
temporal step width t. The stability region separates into an infinite number of disconnected domains. From a practical point of 
view, only the stability domain that includes the point (h,T) - (0,0) is relevant because the temporal step width t is not only 
restricted by the stability criterion (46 1 but also by mc^T/H < n. The later criterion follows from the wave-function's fast phase 
oscillations that are proportional to the typical energy of a wave packet which is about mc^. In stability domains with mc^T/h > n 
the finite difference scheme is stable but does not give accurate numerical solutions. The shape of the stability domains depends 
on the modulus of the vector potential. The larger the vector potential the smaller the stability domains. The region of algebraic 
stability for a finite difference scheme based on five point formulas \26) looks qualitatively similar to the region for three point 
formulas (25 1 but the stability domains are slightly smaller, see Fig. [31 
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FIG. 2: Upper bound r^ax on the temporal step width t as a function of the magnitude of the vector potential Aq. 

It is instructive to analyze the region of algebraic stability A in some limiting cases. Let us consider the case of small step 
widths and small vector potentials with 



n 



<K 1 



mchi 



« 1. 



q\Ao.i\ 



«c 1 



for all / e {1, 



,d} 



first. For second order finite dififerences (25 1 the condition y(f*) < 1 yields in this limit in leading order 



mc^T 1 / inch: \ ^ mc^T 
V /=1 ^ 



< 1 



(47) 



This expression becomes maximal for ^* - +71. For equal grid spacings h\ - . . . - hd - h the condition yiCY' ^ 1 simplifies to 



1 - 2d 



< 1 



which is equivalent to 



A similar calculation for five point finite differences formulas ( [26] l yields the condition 

> V<ic. 

2 T 

The other limit we want to consider is the limit of large spatial step sizes with 

mchi 



(48) 



(49) 



(50) 



n 



» 1 . 



In this limit the condition y{^*) < 1 yields for second order finite differences (25 1 as well as for fourth order finite differences 



( 26 1 the inequality 



mc T 1 mc T I oAq \ mc'T 
cos — — - — sm 



2 h \ mc 



< 1 



(51) 



which no longer depends on f *. The larger the spatial steps the larger the allowed temporal step width. Thus ( |5T| i implies an upper 
bound on the temporal step width which we denote by t^^^. Figure [2] shows Tmax as a function of the magnitude of the vector 
potential Aq. For a vanishing vector potential we have mc^Tmax/^ - ^ and mc^Tmax/fi < ^ otherwise. 



6. Numerical results 



A von-Neumann stability analysis does not account for boundary conditions, it implicitly assumes an infinite grid. Furthermore, 
the von-Neumann stability analysis is not applicable to systems with spatially dependent potentials. Thus, we are going to 
compare our analytical results of Sec.|5]with the stability of actual numerical simulations that work on finite grids and implement 
the boundary condition that the wave function vanishes at border. We will determine the numerical stability for systems with 
spatially dependent potentials as a function of the size of the grid spacings and compare the results with the predictions of our 
von-Neumann stability analysis for homogeneous potentials. In all our numerical simulations atomic units (a.u.) with = 1 a.u., 
q = e - I a.u., m - I a.u., and c ~ 137.036 a.u. will be applied. 
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2nd order finite differences 4th order finite differences 




0.0 0.5 1.0 1.5 2.0 0.0 0.5 1.0 1.5 2.0 0.0 0.5 1.0 1.5 2.0 0.0 0.5 1.0 1.5 2.0 
mch/h mch/ti mch/H mch/h 



FIG. 3: Comparison of the critical lines as predicted by l |46[ l (solid lines) with the points of transition from stability to instability as they have 
been observed in actual numerical simulations (filled circles) of one-dimensional wave packets using finite diff'erences of second order p.5\ and 
of fourth order \26\ , respectively. 



6.1. Systems with homogeneous potentials 



In order to check the condition ( |46| l we determined the region of stability also numerically by propagating a one-dimensional 
wave packet in homogeneous potentials using different grid spacings h and t. The initial wave packet was given by a Gaussian 
superposition of positive energy free particle eigen-states with momentum p 

1 r 1 (1+ ^|TTpy(mcn 1 / (P-Pf , - W], 

^ " V2^ J-.. 2(1 + p2/(mc)2)i/4 \l - VI + py(mc)^) (27rA2)i/4 ^""^ \ 4A^ ^ ' h ) ' ^^^^ 

The initial wave packet has mean momentum p and momentum width A. For our tests we chose p - 2Q a.u., A = 1 a.u. and 
a computational grid that ranges from xiett = -3 a.u. to Xjight — 5 a.u. The wave packet was propagated from time finit - a.u. 
to fend = 1/20 a.u. For a free wave packet not only the pseudo norm (|9 1 is conserved but also the Euclidean norm Thus, 



numerical stability may be checked by determining if the discrete norm (24]i grows exponentially or not as a function of the grid 
spacings h and r. Figure JSlpresents our results for two different vector potentials and finite difference formulas of second order 
( p5| and of fourth order ( |26| , respectively. We observe an excellent correspondence between the von-Neumann stability analysis 
result ([46]l and our numerical findings. 



6.2. Systems with spatially dependent potentials 

We considered so far systems with homogeneous potentials only. In numerical calculations, however, one typically wants 
to simulate systems with spatially dependent potentials to which the technique of the von-Neumann stability analysis is not 
applicable. It is possible to show that the results of von-Neumann stability analysis for systems with homogeneous potentials can 
generally not be extended to inhomogeneous systems. In particular, one can devise finite difference schemes for partial differential 
equations that are instable although the von-Neumann condition ([36]l is satisfied locally everywhere [28 1. 

Assuming that the stability properties of the real space split operator method are well behaved one expects that for systems 
with potentials that vary on length scales that are large compared to the size of the wave packet the stability region is given in 
close approximation by (|46]l. If the stability is actually affected by the potentials' spatial dependence then it is expected to have a 
larger impact, if the potentials vary on the length scale that is given by the wave-function's width. In order to test how a spatially 
dependent vector potential changes the stability properties of the real space split operator method we simulated scattering of a free 
wave packet on a strong laser pulse having a very short wave length. 

We start with a two-dimensional Gaussian wave packet 

w _ 1 r°° r 1 A + vrTp2/(mc)2\ i .xp\ 

- 2nn J_ J_ 2(1 + pVimcyy/^ [l - ^1 + py(mc)^j (InA^y/^ ^""^ \ 4A^ ^ ' n j ' ^^^^ 

The wave packet evolves in a laser pulse of wave length A with the electromagnetic fields 

E(x, t) - Eq sm(k ■ X - cut) fji(k -x - cot), (54a) 
Eo 

B(x, f) = — sin(A; ■ x - tot) fjj{k ■ x - cot) (54b) 
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-0.5 0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 
x/a. u. 

FIG. 4: Left panel: Motion of a wave packet in a ultra strong laser field. Plot shows the center of mass motion (gray line) and a snapshot of the 
density |7]l with the corresponding center of mass (gray circle). Right panel: Stability line for the system in the left panel using fourth order finite 
difi'erences. Dots indicate critical points as determined numerically, the solid line is given by l |46^ using the vector-potential's maximal value. 



with |A:| - In /A - (jjjc and Eq ± k and an envelope function fjj(rf) of total j half cycles including a linear turn-on ramp and a 
linear turn-off ramp of / half cycles and a constant plateau in between, viz. 



'(t] + jn)/il7T) if-j <T]/7r<-j + l, 

1 if + / < 77/7!- < 

-r]/(l7t) if -I < rj/n < 0, 

else. 



(54c) 



The corresponding vector potential follows from ([54| via 



A(x, t) 



f 

Jo 



E{x, f) At' . 



(55) 



For our simulation we chose the parameters A = 1 a. u., /I = 3 a. u., j = i - 2 and \Ao\ - \Eo\/aj - mc/q. The initial wave 
packet has a width of a half Bohr radius and scatters at an ultra strong laser pulse having wave length only six times larger and a 
peak intensity of 1.09 x 10^^ W/cm^. The left panel of Fig.|4]shows the center of mass motion and a snapshot of the density (j?]). 

Determining the numerical stability by propagating the wave function using different step widths h] -h2—h and t we find 
that the stability region for this system with spatially dependent vector potential agrees well with if in ( |46| l the maximal 
value of ( [55] l is used, see right panel of Fig. [4] Numerical instability was detected by determining if \2A) grows exponentially or 
not. Note that for a wave packet in a laser field the norm ( |24| l is not a conserved quantity. Our numerical findings suggest the 
conjecture that the real space split operator method is stable also for inhomogeneous potentials if the von-Neumann condition ( [36] l 
is satisfied locally everywhere. 



7. Conclusions 



We analyzed the stability of a real space split operator method for the time-dependent Klein-Gordon equation and found that 
the method is conditionally stable. For homogeneous potentials the region of algebraic stability was determined as (46 1 via a 
von-Neumann stability analysis. Our numerical simulations support the conjecture that the condition ( |46| is also applicable for 
spatially dependent vector potentials. In the limit of small spatial grid spacings h in all spatial directions and small temporal 
steps T the stability condition becomes Ii/t > Vdc for second order finite differences and ^hlilr) > Vdc for fourth order 
finite differences, respectively. Two conditions that might be more handy for practical applications. Our results may help to 
circumvent numerical instabilities in the application of the real space split operator method for the Klein-Gordon equation by 
choosing appropriate spatial grid spacings and an appropriate temporal step width. 
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